1 C T-00 PROGRAM
2 DIMENSION NDLTH(4),CDH(4),AR(5,5),AI(5,5),BR(5,5),BI(5,5),PN(6),
3 1X(10),Y(10),T(10,10),RE(10,10),IM(10,10),TN(6,10,10)
4 DOUBLE PRECISION THETA,CDH,DLTH,SRMUL,STH,CTH,KR,DKR,PN,FCT,X,Y,BT
5 1,CT,W1,W2,W3,W4,W5,AR,AI,BR,BI,RE,IM,K1,R,DR,ALPHA,A,B,C,D,PI,
6 1THETA1,THETA2,THETA3,THETA4,E,F,G
7 COMPLEX*16 T,TN
8 COMMON/FAC/FCT(100),NRISK,NTRIX
9 NAMELIST/HEJ/AR,AI,BR,BI,RE,IM,T
10 cc CALL UNSTAE
11 PI = 4.0D0*DATAN(1.0D0)
12 NRISK = 57
13 NTRIX = 39
14 NEND = 78
15 NR = 5
16 NP = 0
17 ISYM = 0
18 K1 = 0.5D0
19 A = 1.0D0
20 C = -0.1D0
21 NSECT = 1
22 NDLTH(1) = 192
23 THETA1 = PI
24 CDH(1) = THETA1/DFLOAT(NDLTH(1))
25 N = NRISK
26 L = NTRIX
27 M = NEND-2
28 N1 = N-1
29 DO 5 K = 1,100
30 5 FCT(K) = 0.0D0
31 FCT(1) = 1.0D0
32 DO 6 K = 1,N1
33 6 FCT(K+1) = FCT(K)*DFLOAT(K)
34 FCT(N+1) = 0.1D0**L*FCT(N)*DFLOAT(N)
35 DO 7 K = N,M
36 7 FCT(K+2) = FCT(K+1)*DFLOAT(K+1)
37 ND = 2*NR
38 NR1 = NR+1
39 DO 150 M1 = 1,NR1
40 M = M1-1
41 MD = M
42 IF(M.EQ.0) MD = 1
43 DO 25 J = 1,NR
44 DO 25 I = 1,NR
45 AR(I,J) = 0.0D0
46 AI(I,J) = 0.0D0
47 BR(I,J) = 0.0D0
48 BI(I,J) = 0.0D0
49 25 CONTINUE
50 THETA = 0.0D0
51 DO 90 ISECT = 1,NSECT
52 NTHETA = NDLTH(ISECT)+1
53 DLTH = CDH(ISECT)
54 NFIX = 1
55 DO 89 K = 1,NTHETA
56 IF(K-1)31,31,32
57 31 CONTINUE
58 SRMUL = DLTH*7.0D0/22.5D0
59 IF(ISECT.EQ.1) GO TO 89
60 C Seite 1
61
62
63 GO TO 41
64 32 CONTINUE
65 IF(K-NTHETA)34,33,33
66 33 CONTINUE
67 SRMUL = DLTH*7.0D0/22.5D0
68 IF(ISECT-NSECT)40,89,89
69 34 CONTINUE
70 GO TO (35,36,37,38),NFIX
71 35 CONTINUE
72 SRMUL = DLTH*32.0D0/22.5D0
73 NFIX = 2
74
75
76 GO TO 39
77 36 CONTINUE
78 SRMUL = DLTH*12.0D0/22.5D0
79 NFIX = 3
80 GO TO 39
81 37 CONTINUE
82 SRMUL = DLTH*32.0D0/22.5D0
83 NFIX = 4
84 GO TO 39
85 38 CONTINUE
86 SRMUL = DLTH*14.0D0/22.5D0
87 NFIX = 1
88 39 CONTINUE
89 40 CONTINUE
90 THETA = THETA+DLTH
91 STH = DSIN(THETA)
92 CTH = DCOS(THETA)
93 CALL LEG(THETA,M,NR1,PN)
94 41 CONTINUE
95 CALL TRCIR(STH,CTH,C,A,R,DR)
96 KR = K1*R
97 DKR = K1*DR
98 IF(K.EQ.1) GO TO 52
99 CALL BN(KR,NR1,X,Y)
100 BT= DABS(KR*KR*(X(2)*Y(1)-X(1)*Y(2))-1.0D0)
101 CT= DABS(KR*KR*(X(NR1)*Y(NR)-X(NR)*Y(NR1))-1.0D0)
102 IF(BT-1.D-10)48,48,45
103 45 CONTINUE
104 PRINT 46
105 46 FORMAT ('0 ERROR IN BESSEL TEST')
106 PRINT 47,BT,ISECT,K
107 47 FORMAT(D25.15,2I5)
108 48 CONTINUE
109 IF(CT-1.D-10)52,52,49
110 49 CONTINUE
111 PRINT 50
112 50 FORMAT ('0 ERROR IN NEUMAN TEST')
113 PRINT 51,CT,ISECT,K
114 51 FORMAT(D25.15,2I5)
115 52 CONTINUE
116 DO 88 I = MD,NR
117 DO 88 J = MD,NR
118 I1 = I+1
119 J1 = J+1
120 IF(M.EQ.0) GO TO 74
121 IF(ISYM.GT.0.AND.MOD(I+J,2).EQ.0) GO TO 74
122 C Seite 2
123
124
125 W1 = DFLOAT(I+J)*CTH*PN(I1)*PN(J1)-DFLOAT(I+M)*PN(I)*PN(J1)-DFLOAT
126 1(J+M)*PN(I1)*PN(J)
127 W2 = KR*KR*X(I1)
128
129
130 W3 = W1*W2
131 IF(ISYM.NE.2) GO TO 71
132 IF(J-1)72,71,71
133 71 CONTINUE
134 BI(I,J) = BI(I,J)+SRMUL*W3*Y(J1)/STH
135 72 CONTINUE
136 IF(J-I)74,73,73
137 73 CONTINUE
138 BR(I,J) = BR(I,J)+SRMUL*W3*X(J1)/STH
139 74 CONTINUE
140 IF(ISYM.GT.0.AND.MOD(I+J,2).NE.0) GO TO 84
141 W4 = KR*(KR*X(I)-DFLOAT(I)*X(I1))
142 W5 = DFLOAT(I1)*DKR*STH*X(I1)
143 W1 = PN(I1)*PN(J1)*(W4*(DFLOAT(I*J)*(CTH**2)+DFLOAT(M*M))+DFLOAT
144 1(I*J)*W5*CTH)
145 W2 = DFLOAT(I*(J+M))*PN(I1)*PN(J)*(CTH*W4+W5)
146 W3 = DFLOAT(I*M)*PN(I)*W4*(DFLOAT(J+M)*PN(J)-DFLOAT(J)*CTH*PN(J1))
147 IF(ISYM.NE.2) GO TO 81
148 IF(J-I)82,81,81
149 81 CONTINUE
150 AI(I,J) = AI(I,J)+SRMUL*(W1-W2+W3)*Y(J1)/STH
151
152
153
154 82 CONTINUE
155 IF(J-I)84,83,83
156 83 CONTINUE
157 AR(I,J) = AR(I,J)+SRMUL*(W1-W2+W3)*X(J1)/STH
158 84 CONTINUE
159 88 CONTINUE
160 89 CONTINUE
161 90 CONTINUE
162 DO 92 I = 2,NR
163 JEND = I-1
164 DO 91 J = 1,JEND
165 BR(I,J) = BR(J,I)
166 AR(I,J) = AR(J,I)
167 IF(ISYM.NE.2) GO TO 91
168 BI(I,J) = BI(J,I)
169 AI(I,J) = AI(J,I)
170 91 CONTINUE
171 92 CONTINUE
172 DO 97 I = MD,NR
173 DO 97 J = MD,NR
174 I1 = I+1
175 J1 = J+1
176 W1 = -DSQRT(DFLOAT((2*I+1)*(2*J+1))/DFLOAT(I*I1*J*J1))/2.0D0
177 IF(M)96,96,95
178 95 CONTINUE
179 W1 = W1*DSQRT(FCT(I1-M)*FCT(J1-M)/(FCT(I1+M)*FCT(J1+M)))
180 96 CONTINUE
181 AR(I,J) = W1*AR(I,J)
182 AI(I,J) = W1*AI(I,J)
183 BR(I,J) = DFLOAT(M)*W1*BR(I,J)
184 BI(I,J) = DFLOAT(M)*W1*BI(I,J)
185 97 CONTINUE
186 CALL PERFT(NP,M,NR,ND,AR,AI,BR,BI,T,RE,IM,X,Y)
187 C Seite 3
188
189
190 DO 130 I = 1,ND
191 DO 130 J = 1,ND
192 TN(M1,I,J) = T(I,J)
193 130 CONTINUE
194 IF(M.NE.1) GO TO 150
195 WRITE(6,HEJ)
196
197
198 150 CONTINUE
199 WRITE(11) TN
200 PRINT 161
201 161 FORMAT ('0 TN NOW HOPEFULLY REPLACED INTO DATA SET')
202 STOP
203 cc DEBUG SUBCHK
204 END
205 C Seite 4
206
207
208 C Q-D PROGRAM
209 DIMENSION NDLTH(4),CDH(4),ARR(5,5),AOR(5,5),ARO(5,5),AOO(5,5),
210 1BRR(5,5),BOR(5,5),BRO(5,5),BOO(5,5),CRR(5,5),COR(5,5),CRO(5,5),
211 1COO(5,5),DRR(5,5),DOR(5,5),DRO(5,5),DOO(5,5),PN(7),X1(6),X2(6),
212 1Y1(6),Y2(6),Q(10,10)
213 DOUBLE PRECISION THETA,CDH,DLTH,SRMUL,STH,CTH,PN,DPNI,DPNJ,FCT,B1,
214 1B2,C1,C2,PMY1,PMY2,GA,R,DR,NIND,KAPPA,ALPHA,A,B,C,D,PI,THETA1,
215 1THETA2,THETA3,THETA4,E,F,G
216 COMPLEX*16 CI,K1R,K2R,FMY,DK1R,DK2R,X1,X2,Y1,Y2,Z1I,Z2J,DKRX1I,
217 1DKRY1I,DKRZ1I,DKRX2J,DKRY2J,DKRZ2J,W1,W2,W3,W4,W5,W6,W7,W8,W9,
218 1ARR,AOR,ARO,AOO,BRR,BOR,BRO,BOO,CRR,COR,CRO,COO,DRR,DOR,DRO,DOO,
219 1K1,K2,Q,FC
220 COMMON/FAC/FCT (100),NRISK,NTRIX
221 NAMELIST/HEJ/ARR,AOR,ARO,AOO,BRR,BOR,BRO,BOO,CRR,COR,CRO,COO,DRR,
222 1DOR,DRO,DOO/NEJ/Q
223 cc CALL UNSTAE
224 PI = 4.0D0*DATAN(1.0D0)
225 NRISK = 57
226 NTRIX = 39
227 NEND = 78
228 NR = 5
229 NP = 0
230 LAY = 1
231 PMY1 = 1.0D0
232 PMY2 = 1.0D0
233 K1 = (0.75D0,0.0D0)
234 NIND = 4.0D0/5.0D0
235 KAPPA = 0.0D0
236 C = -0.1D0
237 A = 1.0D0
238 NSECT = 1
239 NDLTH(1) = 192
240 THETA1 = PI
241 CDH(1) = THETA1/DFLOAT(NDLTH(1))
242 N = NRISK
243 L = NTRIX
244 M = NEND-2
245 N1 = N-1
246 DO 5 K = 1,100
247 5 FCT(K) = 0.0D0
248 FCT(1) = 1.0D0
249 DO 6 K = 1,N1
250 6 FCT(K+1) = FCT(K)*DFLOAT(K)
251 FCT(N+1) = 0.1D0**L*FCT(N)*DFLOAT(N)
252 DO 7 K = N,M
253 7 FCT (K+2) = FCT(K+1)*DFLOAT(K+1)
254 FC = (1.0D0,0.0D0)
255 IF(NP.EQ.1) FC = -FC
256 CI = (0.0D0, 1.000)
257 FMY = DCMPLX(PMY1/PMY2,0.0D0)
258 K2 = DCMPLX(NIND,NIND*KAPPA)*K1
259 ND = 2*NR
260 NR1 = NR+1
261 NR2 = NR+2
262 DO 150 M1 = 1,NR1
263 M = M1-1
264 MD = M
265 IF(M.EQ.0) MD = 1
266 DO 25 J = 1,NR
267 DO 25 I = 1,NR
268 C Seite 5
269
270
271 ARR(I,J) = (0.0D0,0.0D0)
272 AOR(I,J) = (0.0D0,0.0D0)
273 ARO(I,J) = (0.0D0,0.0D0)
274 AOO(I,J) = (0.0D0,0.0D0)
275 BRR(I,J) = (0.0D0,0.0D0)
276 BOR(I,J) = (0.0D0,0.0D0)
277 BRO(I,J) = (0.0D0,0.0D0)
278 BOO(I,J) = (0.0D0,0.0D0)
279 CRR(I,J) = (0.0D0,0.0D0)
280 COR(I,J) = (0.0D0,0.0D0)
281 CRO(I,J) = (0.0D0,0.0D0)
282 COO(I,J) = (0.0D0,0.0D0)
283 DRR(I,J) = (0.0D0,0.0D0)
284 DOR(I,J) = (0.0D0,0.0D0)
285 DRO(I,J) = (0.0D0,0.0D0)
286 DOO(I,J) = (0.0D0,0.0D0)
287 25 CONTINUE
288 THETA = 0.0D0
289 DO 90 ISECT = 1,NSECT
290 NTHETA = NDLTH(ISECT)+1
291 DLTH = CDH(ISECT)
292 NFIX = 1
293 DO 89 K = 1,NTHETA
294 IF(K-1)31,31,32
295 31 CONTINUE
296 SRMUL = DLTH*7.0D0/22.5D0
297 IF(ISECT.EQ.1) GO TO 89
298 GO TO 41
299 32 CONTINUE
300 IF(K-NTHETA)34,33,33
301 33 CONTINUE
302 SRMUL = DLTH*7.0D0/22.5D0
303 IF(ISECT-NSECT)40,89,89
304 34 CONTINUE
305 GO TO (35,36,37,38),NFIX
306 35 CONTINUE
307 SRMUL = DLTH*32.0D0/22.5D0
308 NFIX = 2
309 GO TO 39
310 36 CONTINUE
311 SRMUL = DLTH*12.0D0/22.5D0
312 NFIX = 3
313 GO TO 39
314 37 CONTINUE
315 SRMUL = DLTH*32.0D0/22.5D0
316 NFIX = 4
317 GO TO 39
318 38 CONTINUE
319 SRMUL = DLTH*14.0D0/22.5D0
320 NFIX = 1
321 39 CONTINUE
322 40 CONTINUE
323 THETA = THETA+DLTH
324 STH = DSIN(THETA)
325 CTH = DCOS(THETA)
326 CALL LEG(THETA,M,NR2,PN)
327 C Seite 6
328
329
330 41 CONTINUE
331 CALL TRCIR(STH,CTH,C,A,R,DR)
332 K1R= K1*DCMPLX(R,0.0D0)
333 DK1R = K1*DCMPLX(DR,0.0D0)
334
335
336 K2R = K2*DCMPLX(R,0.0D0)
337 IF(K.EQ.1) GO TO 52
338 CALL CBN(K1R,NR1,X1,Y1)
339 CALL CBN(K2R,NR1,X2,Y2)
340 B1 = CDABS(K1R*K1R*(X1(2)*Y1(1)-X1(1)*Y1(2))-DCMPLX(1.0D0,0.0D0))
341 B2 = CDABS(K2R*K2R*(X2(2)*Y2(1)-X2(1)*Y2(2))-DCMPLX(1.0D0,0.0D0))
342 C1 = CDABS(K1R*K1R*(X1(NR1)*Y1(NR)-X1(NR)*Y1(NR1))-DCMPLX(1.0D0,0.
343 10D0))
344 C2 = CDABS(K2R*K2R*(X2(NR1)*Y2(NR)-X2(NR)*Y2(NR1))-DCMPLX(1.0D0,0.
345 10D0))
346 IF(B1.LT.1.D-10.AND.C1.LT.1.D-10) GO TO 47
347 PRINT 45
348 45 FORMAT('0 ERROR IN BESSEL-NEUMAN-1 TEST')
349 PRINT 46, B1,C1,ISECT,K
350 46 FORMAT(2D25.15,2I5)
351 47 CONTINUE
352 IF(B2.LT.1.D-10.AND.C2.LT.1.D-10) GO TO 52
353 PRINT 48
354 48 FORMAT ('0 ERROR IN BESSEL-NEUMAN-2 TEST')
355 49 FORMAT(2D25.15,2I5)
356 PRINT 49,B2,C2,ISECT,K
357 52 CONTINUE
358 DO 88 I = MD,NR
359 DO 88 J = MD,NR
360 I1 = I+J
361 J1 = J+1
362 I2 = I + 2
363 J2 = J+2
364 DPNI = -(DFLOAT(I1)*CTH*PN(I1)-DFLOAT(I1-M)*PN(I2))/STH
365 DPNJ = -(DFLOAT(J1)*CTH*PN(J1)-DFLOAT(J1-M)*PN(J2))/STH
366 Z1I = X1(I1)+CI*Y1(I1)
367 DKRX1I = K1R*X1(I)-DCMPLX(DFLOAT(I),0.0D0)*X1(I1)
368 DKRY1I = K1R*Y1(I)-DCMPLX(DFLOAT(I),0.0D0)*Y1(I1)
369 DKRZ1I = DKRX1I+CI*DKRY1I
370 Z2J = X2(J1)*CI*Y2(J1)
371 DKRX2J = K2R*X2(J)-DCMPLX(DFLOAT(J),0.0D0)*X2(J1)
372 DKRY2J = K2R*Y2(J)-DCMPLX(DFLOAT(J),0.0D0)*Y2(J1)
373 DKRZ2J = DKRX2J+CI*DKRY2J
374 W1 = DCMPLX(SRMUL*STH*(DPNI*DPNJ*DFLOAT(M*M)*PN(I1)*PN(J1)/STH**2)
375 1,0.0D0)
376 W2 = DCMPLX(SRMUL*STH*DFLOAT(I*I1)*PN(I1)*DPNJ,0.0D0)
377 W3 = DCMPLX(SRMUL*STH*DFLOAT(J*J1)*PN(J1)*DPNI,0.0D0)
378 ARR(I,J) = ARR(I,J)+K1R*DKRX1I*X2(J1)*W1+DK1R*X1(I1)*X2(J1)*W2
379 AOR(I,J) = AOR(I,J)+K1R*DKRZ1I*X2(J1)*W1+DK1R*Z1I*X2(J1)*W2
380 DRR(I,J) = DRR(I,J)+K1R*X1(I1)*DKRX2J*W1+DK1R*X1(I1)*X2(J1)*W3
381 DOR(I,J) = DOR(I,J)+K1R*Z1I*DKRX2J*W1+DK1R*Z1I*X2(J1)*W3
382 IF(LAY.EQ.0) GO TO 53
383 ARO(I,J) = ARO(I,J)+K1R*DKRX1I*Z2J*W1+DK1R*X1(I1)*Z2J*W2
384 AOO(I,J) = AOO(I,J)+K1R*DKRZ1I*Z2J*W1*DK1R*Z1I*Z2J*W2
385 DRO(I,J) = DRO(I,J)+K1R*X1(I1)*DKRZ2J*W1+DK1R*X1(I1)*Z2J*W3
386 DOO(I,J) = DOO(I,J)+K1R*Z1I*DKRZ2J*W1+DK1R*Z1I*Z2J*W3
387 53 CONTINUE
388 IF(M.EQ.0) GO TO 88
389 C Seite 7
390
391
392 W4 = DCMPLX(SRMUL*(DPNI*PN(J1)+PN(I1)*DPNJ),0.0D0)
393 W5 = DCMPLX(SRMUL*PN(I1)*PN(J1),0.0D0)
394 BRR(I,J) = BRR(I,J)+K1R*K1R*X1(I1)*X2(J1)*W4
395 BOR(I,J) = BOR(I,J)+K1R*K1R*Z1I*X2(J1)*W4
396 W6 = DK1R*(X1(I1)*DKRX2J*DCMPLX(DFLOAT(I*(I+1)),0.0D0)+DCMPLX(
397 1DFLOAT(J*(J+1)),0.0D0)*DKRX1I*X2(J1))/K1R
398 CRR(I,J) = CRR(I,J)+DKRX1I*DKRX2J*W4+W5*W6
399 W6 = DK1R*(Z1I*DKRX2J*DCMPLX(DFLOAT (I*(I+1)),0.0D0)+DCMPLX(
400 1DFLOAT(J*(J+1)),0.0D0)*DKRZ1I*X2(J1))/K1R
401 COR(I,J) = COR(I,J)+DKRZ1I*DKRX2J*W4+W5*W6
402 IF(LAY.EQ.0) GO TO 88
403 BRO(I,J) = BRO(I,J)+K1R*K1R*X1(I1)*Z2J*W4
404 BOO(I,J) = BOO(I,J)+K1R*K1R*Z1I*Z2J*W4
405 W6 = DK1R*(X1(I1)*DKRZ2J*DCMPLX(DFLOAT(I*(I+1)),0.0D0)+DCMPLX(
406 1DFLOAT(J*(J+1)),0.0D0)*DKRX1I*Z2J)/K1R
407 CRO(I,J) = CRO(I,J)+DKRX1I*DKRZ2J*W4+W5*W6
408 W6 = DK1R*(Z1I*DKRZ2J*DCMPLX(DFLOAT(I*(I+1)),0.0D0)+DCMPLX(
409 1DFLOAT(J*(J+1)),0.0D0)*DKRZ1I*Z2J)/K1R
410 COO(I,J) = COO(I,J)+DKRZ1I*DKRZ2J*W4+W5*W6
411 88 CONTINUE
412 89 CONTINUE
413 90 CONTINUE
414 DO 98 I = MD,NR
415 DO 98 J = MD,NR
416 I1 = I+1
417 J1 = J+1
418 GA = DSQRT(DFLOAT((2*I+1)*(2*J+1))/DFLOAT(I*I1*J*J1))/2.0D0
419 IF(M)96,96,95
420 95 CONTINUE
421 GA = GA*DSQRT(FCT(I1-M)*FCT(J1-M)/FCT(I1+M)*FCT(J1*M))
422 96 CONTINUE
423 W1 = DCMPLX(GA,0.0D0)
424 ARR(I,J) = W1*ARR(I,J)
425 AOR(I,J) = W1*AOR(I,J)
426 DRR(I,J) = W1*DRR(I,J)
427 DOR(I,J) = W1*DOR(I,J)
428 IF(LAY.EQ.0) GO TO 97
429 ARO(I,J) = W1*ARO(I,J)
430 AOO(I,J) = W1*AOO(I,J)
431 DRO(I,J) = W1*DRO(I,J)
432 DOO(I,J) = W1*DOO(I,J)
433 97 CONTINUE
434 IF(M.EQ.0) GO TO 98
435 W1 = DCMPLX(DFLOAT(M),0.0D0)*W1
436 BRR(I,J) = W1*BRR(I,J)
437 BOR(I,J) = W1*BOR(I,J)
438 CRR(I,J) = W1*CRR(I,J)
439 COR(I,J) = W1*COR(I,J)
440 IF(LAY.EQ.0) GO TO 98
441 BRO(I,J) = W1*BRO(I,J)
442 BOO(I,J) = W1*BOO(I,J)
443 CRO(I,J) = W1*CRO(I,J)
444 COO(I,J) = W1*COO(I,J)
445 98 CONTINUE
446 IF(M.GT.3) GO TO 99
447 WRITE(6,HEJ)
448 99 CONTINUE
449 C Seite 8
450
451
452 LLAY = 2+LAY*2
453 DO 120 LL = 1,LLAY
454 DO 110 I = 1,NR
455 DO 110 J = 1,NR
456 GO TO (101,102,103,104), LL
457 101 CONTINUE
458 Q(2*I-1,2*J-1) = -ARR(I,J)+FMY*DRR(I,J)
459 Q(2*I-1,2*J) = FC*(K1*CRR(I,J)/K2+K2*FMY*BRR(I,J)/K1)
460 Q(2*I,2*J-1) = -FC*(BRR(I,J)+FMY*CRR(I,J))
461 Q(2*I,2*J) = (K1*DRR(I,J)/K2-K2*FMY*ARR(I,J)/K1)
462
463
464 GO TO 110
465 102 CONTINUE
466 Q(2*I-1,2*J-1) =-AOR(I,J)+FMY*DOR(I,J)
467 Q(2*I-1,2*J) = FC*K1*COR(I,J)/K2+K2*FMY*BOR(I,J)/K1
468 Q(2*I,2*J-1) = -FC*(BOR(I,J)*FMY*COR(I,J))
469 Q(2*I,2*J) = (K1*DOR(I,J)/K2-K2*FMY*AOR(I,J)/K1)
470 GO TO 110
471 103 CONTINUE
472 Q(2*I-1,2*J-1) = -ARO(I,J)+FMY*DRO(I,J)
473 Q(2*I-1,2*J) = FC*(K1*CRO(I,J)/K2+K2*FMY*BRO(I,J)/K1)
474 Q(2*I,2*J-1) = -FC*(BRO(I,J)+FMY*CRO(I,J))
475 Q(2*I,2*J) = (K1*DRO(I,J)/K2-K2*FMY*AOO(I,J)/K1)
476 104 CONTINUE
477 Q(2*I-1,2*J-1) = -AOO(I,J)+FMY*DOO(I,J)
478 Q(2*I-1,2*J) = FC*(K1*COO(I,J)/K2+K2*FMY*BOO(I,J)/K1)
479 Q(2*I,2*J-1) = -FC*(BOO(I,J)+FMY*COO(I,J))
480 Q(2*I,2*J) = (K1*DOO(I,J)/K2-K2*FMY*AOO(I,J)/K1)
481 110 CONTINUE
482 WRITE (6,NEJ)
483 GO TO (111,111,112,112)LL
484 111 CONTINUE
485 WRITE (21) Q
486 END FILE 21
487 GO TO 120
488 112 CONTINUE
489 WRITE(22) Q
490 END FILE 22
491 120 CONTINUE
492 PRINT 777, M1
493 777 FORMAT(I4)
494 150 CONTINUE
495 PRINT 199
496 199 FORMAT('0 Q-MATRICES NOW HOPEFULLY PLACED INTO TAPE')
497 STOP
498 cc DEBUG SUBCHK
499 END
500 C Seite 9
501
502
503 C Q-T PROGRAM
504 DIMENSION QRR(10,10),QOR(10,10),LL(10),MM(10),TN(6,10,10)
505 COMPLEX*16 QRR,QOR,TN,P,D
506 cc CALL UNSTAE
507 NR = 5
508 NR1 = NR*1
509 ND = 2*NR
510 DO 30 M1 = 1,NR1
511 READ(21,END = 8) QRR
512 GO TO 9
513 8 READ(21,END=999) QRR
514 9 CONTINUE
515 READ(21,END=12) QOR
516 GO TO 13
517 12 READ(21,END=999) QOR
518 13 CONTINUE
519 IF(M1.LT.3) GO TO 16
520 MD1 = 2*M1-4
521 DO 15 I = 1,MD1
522 QOR = (1.0D0,0.0D0)
523 15 CONTINUE
524 16 CONTINUE
525 CALL MCNV(QOR,ND,D,LL,MM)
526 DO 21 I = 1,ND
527 DO 21 J = 1,ND
528 P = (0.0D0,0.0D0)
529 DO 20 K = 1,ND
530 P = P-QRR(I,K)*QOR(K,J)
531 20 CONTINUE
532 TN(M1,I,J) = P
533 21 CONTINUE
534 30 CONTINUE
535 WRITE (11) TN
536 PRINT 50
537 50 FORMAT('0 TN NOW HOPEFULLY REPLACED INTO DATASE')
538 999 CONTINUE
539 STOP
540 cc DEBUG SUBCHK
541 END
542 C Seite 10
543
544
545 C Q, T-T PROGRAM
546 DIMENSION Q1 (10,10),Q2(10,10),Q3(10,10),T(10,10),LL(10),MM(10),TN
547 1(6,10,10)
548 COMPLEX*16 Q1,Q2,Q3,T,P,D,TN
549 cc CALL UNSTAE
550 NR = 5
551 NR1 = NR+1
552 ND = 2*NR
553 DO 50 M1 = 1,NR1
554 READ(21,END=8) Q1
555 GO TO 9
556 8 READ(21,END=999) Q1
557 9 CONTINUE
558 READ(21,END=12) Q2
559 GO TO 13
560 12 READ(21,END=999) Q2
561 13 CONTINUE
562 READ(22,END=16) Q3
563 GO TO 17
564 16 READ(22,END=999) Q3
565 17 CONTINUE
566 READ(23,END=20) T
567 GO TO 21
568 20 READ(23,END=999) T
569 21 CONTINUE
570 DO 25 J = 1,ND
571 DO 25 I = 1,ND
572 P = (0.0D0,0.0D0)
573 DO 24 K = 1,ND
574 P = P+Q3(I,K)*T(K,J)
575 24 CONTINUE
576 Q1(I,J) = Q1(I,J)+P
577 25 CONTINUE
578 READ(22,END=30) Q3
579 GO TO 31
580 30 READ(22,END=999) Q3
581 31 CONTINUE
582 DO 36 J = 1,ND
583 DO 36 I = 1,ND
584 P = (0.0D0,0.0D0)
585 DO 35 K = 1,ND
586 P = P+Q3(I,K)*T(K,J)
587 35 CONTINUE
588 Q2(I,J) = Q2(I,J)+P
589 36 CONTINUE
590 IF(M1.LT.3) GO TO 41
591 MD1 = 2*M1-4
592 DO 40 I = 1,MD1
593 Q2(I,I) = (1.0D0,0.0D0)
594 40 CONTINUE
595 41 CONTINUE
596 CALL MCNV(Q2,ND,D,LL,MM)
597 DO 46 J = 1,ND
598 DO 46 I = 1,ND
599 P = (0.0D0,0.0D0)
600 DO 45 K = 1,ND
601 P = P+Q1(I,K)*Q2(K,J)
602 45 CONTINUE
603 T(I,J) = -P
604 TN(M1,I,J) = -P
605 C Seite 11
606
607
608 46 CONTINUE
609 50 CONTINUE
610 WRITE(11) TN
611 PRINT 70
612 70 FORMAT('0 TN NOW HOPEFULLY REPLACED INTO DATASET')
613 999 CONTINUE
614 STOP
615 cc DEBUG SUBCHK
616 END
617 C Seite 12
618
619
620 C TN-T PROGRAM FOR Q,T-T PROGRAM
621 DIMENSION TN(6,10,10),T(10,10)
622 COMPLEX*16 TN,T
623 cc CALL UNSTAE
624 ND = 10
625 NS = ND/2+1
626 READ(11) TN
627 DO 20 M = 1,NS
628 DO 10 I = 1,ND
629 DO 10 J = 1,ND
630 T(I,J) = TN(M,I,J)
631 10 CONTINUE
632 WRITE(21) T
633 END FILE 21
634 20 CONTINUE
635 STOP
636 cc DEBUG SUBCHK
637 END
638 C Seite 13
639
640
641 C T1,T2-T PROGRAM
642 DIMENSION R1(10,10),R2(10,10),R3(10,10),R4(10,10),AT1(10,10),AT2(
643 110,10),BT1(10,10),BT2(10,10),T(10,10),RET(10,10),RETT(10,10),TRAN(
644 110,10),T1(10,10),T2(10,10),X1(11),X2(11),Y1(11),Y2(11),LL(10),MM(
645 110),TN(6,10,10)
646 COMPLEX*16 R1,R2,R3,R4,AT1,AT2,BT1,BT2,T,RET,RETT,TRAN,TN,T1,T2,
647 1W1,W2,W3,W4,S,DD
648 DOUBLE PRECISION X1,X2,Y1,Y2,SEP,FCT
649 COMMON/FAC/FCT(100),NRISK,NTRIX
650 NAMELIST/HEJ/T,RET,RETT,TRAN
651 CALL UNSTAE
652 NRISK = 57
653 NTRIX = 39
654 NEND = 78
655 ND = 10
656 NP = 0
657 SEP = 1.5D0
658 N = NRISK
659 L = NTRIX
660 M = NEND-2
661 N1 = N-1
662 DO 5 K = 1,100
663 5 FCT(K) = 0.0D0
664 FCT(1) = 1.0D0
665 DO 6 K=1,N1
666 6 FCT(K+1) = FCT(K)*DFLOAT(K)
667 FCT(N+1) = 0.1D0**L*FCT(N)*DFLOAT(N)
668 DO 7 K = N,M
669 7 FCT(K+2) = FCT(K+1)*DFLOAT(K+1)
670 NS = ND/2+1
671 DO 200 M1 = 1,NS
672 M = M1-1
673 READ(21,END=8) T1
674 GO TO 9
675 8 READ(21,END=999) T1
676 9 CONTINUE
677 READ(22,END=12) T2
678 GO TO 13
679 12 READ(22,END=999)
680 13 CONTINUE
681 IF(M-1)20,20,21
682 20 CONTINUE
683 MD = 1
684 GO TO 22
685 21 CONTINUE
686 MD = 2*M-1
687 22 CONTINUE
688 DO 34 I = 1,ND
689 DO 34 J = 1,ND
690 T(I,J) = (0.0D0,0.0D0)
691 TN(M1,I,J) = (0.0D0,0.0D0)
692 IF(I-J) 30,31,30
693 30 CONTINUE
694 AT1(I,J) = (0.0D0,0.0D0)
695 AT2(I,J) = (0.0D0,0.0D0)
696 BT1(I,J) = (0.0D0,0.0D0)
697 BT2(I,J) = (0.0D0,0.0D0)
698 GO TO 34
699 31 CONTINUE
700 C Seite 14
701
702 IF(M-2)30,32,32
703 32 CONTINUE
704
705
706 MD1 = MD-1
707 IF(MD1-J)30,33,33
708 33 CONTINUE
709 AT1(I,J) = (1.0D0,0.0D0)
710 AT2(I,J) = (1.0D0,0.0D0)
711 BT1(I,J) = (1.0D0,0.0D0)
712 BT2(I,J) = (1.0D0,0.0D0)
713 34 CONTINUE
714 CALL VR(SEP,NP,M,ND,R1,R2,R3,R4,X1,X2,Y1,Y2)
715 DO 41 I = MD,ND
716 DO 41 J = MD,ND
717 W1 = (0.0D0,0.0D0)
718 W2 = (0.0D0,0.0D0)
719 W3 = (0.0D0,0.0D0)
720 W4 = (0.0D0,0.0D0)
721 DO 40 K = MD,ND
722 W1 = W1+R2(I,K)*T2(K,J)
723 W2 = W2+R3(I,K)*T1(K,J)
724 W3 = W3+R4(I,K)*T1(K,J)
725 W4 = W4+R4(I,K)*T2(K,J)
726 40 CONTINUE
727 T(I,J) = W1
728 RET(I,J) = W2
729 RETT(I,J) = W3
730 TRAN(I,J) = W4
731 41 CONTINUE
732 DO 45 I = MD,ND
733 DO 45 J = MD,ND
734 W1 = (0.0D0,0.0D0)
735 W2 = (0.0D0,0.0D0)
736 W3 = (0.0D0,0.0D0)
737 W4 = (0.0D0,0.0D0)
738 DO 42 K = MD,ND
739 W1 = W1+T(I,K)*RET(K,J)
740 W2 = W2+RET(I,K)*T(K,J)
741 W3 = W3+T(I,K)*R1(J,K)
742 W4 = W4+RET(I,K)*R1(K,J)
743 42 CONTINUE
744 IF(I-J)44,43,44
745 43 CONTINUE
746 AT1(I,J) = (1.0D0,0.0D0)-W1
747 AT2(I,J) = (1.0D0,0.0D0)-W2
748 BT1(I,J) = (1.0D0,0.0D0)+W3
749 BT2(I,J) = (1.0D0,0.0D0)+W4
750 GO TO 45
751 44 CONTINUE
752 AT1(I,J) = -W1
753 AT2(I,J) = -W2
754 BT1(I,J) = W3
755 BT2(I,J) = W4
756 45 CONTINUE
757 CALL MCNV(AT1,ND,DD,LL,MM)
758 CALL MCNV(AT2,ND,DD,LL,MM)
759 DO 47 I = MD,ND
760 DO 47 J = MD,ND
761 C Seite 15
762
763
764 W1 = (0.0D0,0.0D0)
765 W2 = (0.0D0,0.0D0)
766 W3 = (0.0D0,0.0D0)
767 W4 = (0.0D0,0.0D0)
768 DO 46 K = MD,ND
769
770
771 W1 = W1+RETT(I,K)*AT1(K,J)
772 W2 = W2+TRAN(I,K)*AT2(K,J)
773 W3 = W3+BT1(I,K)*R4(K,J)
774 W4 = W4+BT2(I,K)*R4(J,K)
775 46 CONTINUE
776 R1(I,J) = W1
777 R2(I,J) = W2
778 R3(I,J) = W3
779 RET(I,J) = W4
780 47 CONTINUE
781 DO 70 I = MD,ND
782 DO 70 J = MD,ND
783 W1 = (0.0D0,0.0D0)
784 W2 = (0.0D0,0.0D0)
785 DO 49 K = MD,ND
786 W1 = W1+R1(I,K)*R3(K,J)
787 W2 = W2+R2(I,K)*RET(K,J)
788 49 CONTINUE
789 S = W1+W2
790 T(I,J) = S
791 TN(M1,I,J) = S
792 70 CONTINUE
793 IF((M-1).NE.0) GO TO 200
794 DO 81 I = 1,ND
795 DO 81 J = 1,ND
796 S = (0.0D0,0.0D0)
797 DO 80 K = MD,ND
798 S = S-DCONJG(T(K,I))*T(K,J)
799 80 CONTINUE
800 RET(I,J) = S
801 RETT(I,J) = T(I,J)-T(J,I)
802 TRAN(I,J)=T(I,J)-T(J,I)
803 81 CONTINUE
804 WRITE(6,HEJ)
805 200 CONTINUE
806 WRITE(11) TN
807 PRINT 98
808 98 FORMAT('0 TN NOW HOPEFULLY REPLACED INTO DATASET')
809 999 CONTINUE
810 STOP
811 cc DEBUG SUBCHK
812 END
813 C Seite 16
814
815
816
817 C TN-T PROGRAM NR 1 FOR T1,T2-T PROGRAM
818 DIMENSION TN(6,10,10),T(10,10)
819 COMPLEX*16 TN,T
820 cc CALL UNSTAE
821 ND = 10
822 NS = ND/2+1
823 READ(11) TN
824 DO 20 M = 1,NS
825 DO 10 I = 1,ND
826 DO 10 J = 1,ND
827 T(I,J) = TN(M,I,J)
828 10 CONTINUE
829 WRITE(21) T
830 END FILE 21
831 20 CONTINUE
832 STOP
833 cc DEBUG SUBCHK
834 END
835 C Seite 17
836
837
838 C TN-T PROGRAM NR 2 FOR T1,T2-T PROGRAM
839 DIMENSION TN(6,10,10),T(10,10)
840 COMPLEX*16 TN,T
841 cc CALL UNSTAE
842 ND = 10
843 NS = ND/2+1
844 READ(11) TN
845 DO 20 M = 1,NS
846 DO 10 I = 1,ND
847 DO 10 J = 1,ND
848 T(I,J) = TN(M,I,J)
849 10 CONTINUE
850 WRITE(21) T
851 END FILE 21
852 20 CONTINUE
853 STOP
854 cc DEBUG SUBCHK
855 END
856 C Seite 18
857
858
859 C T-T PROGRAM
860 DIMENSION TN(6,10,10),T(10,10),R1(10,10),R2(10,10),R3(10,10),R4
861 1(10,10),X1(11),X2(11),Y1(11),Y2(11)
862 DOUBLE PRECISION TR,X1,X2,Y1,Y2,FCT
863 COMPLEX*16 TN,T,R1,R2,R3,R4,P
864 COMMON/FAC/FCT(100),NRISK,NTRIX
865 NAMELIST/HEJ/T
866 NRISK = 57
867 NTRIX = 39
868 NEND = 78
869 ND = 10
870 NP = 0
871 TR = 0.5D0
872 N = NRISK
873 L = NTRIX
874 M = NEND-2
875 N1 = N-1
876 DO 5 K = 1,100
877 5 FCT(K) = 0.0D0
878 FCT(1) = 1.0D0
879 DO 6 K = 1,N1
880 6 FCT(K+1) = FCT(K)*DFLOAT(K)
881 FCT(N+1) = 0.1D0**L*FCT(N)*DFLOAT(N)
882 DO 7 K = N,M
883 7 FCT(K+2) = FCT(K+1)*DFLOAT(K+1)
884 READ(11) TN
885 NS = ND/2+1
886 DO 20 M1 = 1,NS
887 M = M1-1
888 CALL VR(TR,NP,M,ND,R1,R2,R3,R4,X1,X2,Y1,Y2)
889 DO 11 J = 1,ND
890 DO 11 I = 1,ND
891 P = (0.0D0,0.0D0)
892 DO 10 K = 1,ND
893 P = P+TN(M1,I,K)*R1(J,K)
894 10 CONTINUE
895 R4(I,J) = P
896 11 CONTINUE
897 DO 16 J = 1,ND
898 DO 16 I = 1,ND
899 P = (0.0D0,0.0D0)
900 DO 15 K = 1,ND
901 P = P+R1(I,K)*R4(K,J)
902 15 CONTINUE
903 T(I,J) = P
904 TN(M1,I,J) = P
905 16 CONTINUE
906 WRITE(6,HEJ)
907 20 CONTINUE
908 WRITE(12) TN
909 PRINT 70
910 70 FORMAT('0 TN NOW HOPEFULLY REPLACED INTO DATSET')
911 STOP
912 cc DEBUG SUBCHK
913 END
914 C Seite 19
915
916
917 C PSIS PROGRAM
918 DIMENSION TN(6,10,10),AP(10),PN(7),FEXP(10),PSIR(10),PSITH(10),
919 1PSIFH(10),X(7),Y(7)
920 DOUBLE PRECISION BHETA,PN,DIST,THETA,FHI,X,Y,FCT,DIFSC,RAD,AMPL,
921 1BHETAD,FHID,THETAD,PI,KV,AMPL2,A,B,C,D,E
922 COMPLEX*16 TN,AP,PSIR,PSITH,PSIFH,FEXP,Q1,Q2,P1,P2,R1,R2,FC
923 COMMON/FAC/FCT(100),NRISK,NTRIX
924 NAMELIST/HEJ/Q1,Q2,P1,P2,R1,R2
925 CALL UNSTAE
926 PI = 4.0D0*DATAN(1.0D0)
927 NRISK = 57
928 NTRIX = 39
929 NEND = 78
930 NP = 1
931 NPCHAN = 1
932 ND = 10
933 NB = 1
934 KV = 2.0D0*PI/15.53333D0
935 BHETA = PI/2.0D0
936 BHETAD = BHETA*180.0D0/PI
937 DIST = 1.0D0
938 THETA = PI/2.0D0
939 THETAD = THETA*180.0D0/PI
940 FHI = 0.0D0
941 FHID = FHI*180.0D0/PI
942 N = NRISK
943 L = NTRIX
944 M = NEND-2
945 N1 = N-1
946 DO 5 K = 1,100
947 5 FCT(K) = 0.0D0
948 FCT(1) = 1.0D0
949 DO 6 K = 1,N1
950 6 FCT(K+1) = FCT(K)*DFLOAT(K)
951 FCT(N+1) = 0.1D0**L*FCT(N)*DFLOAT(N)
952 DO 7 K = N,M
953 7 FCT (K+2) = FCT(K+1)*DFLOAT(K+1)
954 NPC = NPCHAN+1
955 GO TO (11,12), NPC
956 11 CONTINUE
957 FC = (1.0D0,0.0D0)
958 GO TO 13
959 12 CONTINUE
960 FC = (-1.0D0,0.0D0)
961 13 CONTINUE
962 READ (11) TN
963 NS = ND/2+1
964 P1 = (0.0D0,0.0D0)
965 P2 = (0.0D0,0.0D0)
966 DO 23 M1 = 1,NS
967 M = M1-1
968 MD = 2*M-1
969 IF((M-2).LT.0) MD = 1
970 CALL VKOEF(BHETA,NP,M,ND,AP,PN)
971 CALL VPSI(DIST,THETA,FHI,NP,NB,M,ND,PSIR,PSITH,PSIFH,PN,X,Y)
972 DO 21 I = MD,ND
973 Q1 = (0.0D0,0.0D0)
974 DO 20 J = MD,ND
975 Q1 = Q1+TN(M1,I,J)*AP(J)*FC**(I+J)
976 C Seite 20
977
978
979 20 CONTINUE
980
981
982 FEXP(I) = Q1
983 21 CONTINUE
984 Q1 = (0.0D0,0.0D0)
985 Q2 = (0.0D0,0.0D0)
986 DO 22 I = MD,ND
987 Q1 = Q1+PSITH(I)*FEXP(I)
988 Q2 = Q2+PSIFH(I)*FEXP(I)
989 22 CONTINUE
990 P1 = P1+Q1
991 P2 = P2+Q2
992 R1 = P1/DCMPLX(KV,0.0D0)
993 R2 = P2/DCMPLX(KV,0.0D0)
994 WRITE(6,HEJ)
995 AMPL2= CDABS(P1)**2+CDABS(P2)**2
996 AMPL = DSQRT(AMPL2)
997 A = THETAD
998 B = FHID
999 C = AMPL
1000 D = AMPL2
1001 E = AMPL/KV
1002 PRINT 40,A,B,C,D,E
1003 40 FORMAT(5D25.15)
1004 23 CONTINUE
1005 100 CONTINUE
1006 300 CONTINUE
1007 999 CONTINUE
1008 STOP
1009 cc DEBUG SUBCHK
1010 END
1011 C Seite 21
1012
1013
1014 C T TEST PROGRAM
1015 DIMENSION TN(6,20,20),T(10,10),RET(10,10),RETT(10,10),TRAN(10,10)
1016 COMPLEX*16 TN,T,RET,RETT,TRAN,S
1017 NAMELIST/HEJ/T,RET,RETT,TRAN
1018 CALL UNSTAE
1019 ND = 10
1020 READ (11) TN
1021 NS = ND/2+1
1022 DO 20 M = 1,NS
1023 DO 10 I = 1,ND
1024 DO 10 J = 1,ND
1025 T(I,J) = TN(M,I,J)
1026 10 CONTINUE
1027 DO 12 I = 1,ND
1028 DO 12 J = 1,ND
1029 S = (0.0D0,0.0D0)
1030 DO 11 K = 1,ND
1031 S = S-DCONJG(T(K,I))*T(K,J)
1032 11 CONTINUE
1033 RET(I,J) = S
1034 RETT(I,J) = T(I,J)-S
1035 TRAN(I,J) = T(I,J)-T(J,I)
1036 12 CONTINUE
1037 WRITE(6,HEJ)
1038 IF(M.GT.3) GO TO 99
1039 20 CONTINUE
1040 99 CONTINUE
1041 STOP
1042 cc DEBUG SUBCHK
1043 END
1044 C Seite 22
1045
1046
1047 SUBROUTINE BESSEL(NORDER,ARGMNT,ANSWR,IERROR)
1048 DOUBLE PRECISION ARGMNT,ANSWR,X,SUM,APR,TOPR,CI,CNI,ACR,PROD,FACT
1049 IERROR = 0
1050 N = NORDER
1051 X = ARGMNT
1052 SUM = 1.0D0
1053 APR = 1.0D0
1054 TOPR = -0.5D0*X*X
1055 CI = 1.0D0
1056 CNI = DFLOAT(2*N+3)
1057 DO 60 I = 1,100
1058 ACR = TOPR*APR/(CI*CNI)
1059 SUM = SUM+ACR
1060 IF(DABS(ACR/SUM)-1.0D-20)100,100,40
1061 40 APR = ACR
1062 CI = CI+1.0D0
1063 CNI = CNI+2.0D0
1064 60 CONTINUE
1065 IERROR = I
1066 PRINT 10
1067 10 FORMAT(24HO ERROR IN SUM OF BESSEL)
1068 GO TO 200
1069 100 PROD = DFLOAT(2*N+1)
1070 FACT = 1.0D0
1071 IF(N)160,160,120
1072 120 DO 140 IFCT = 1,N
1073 FACT = FACT*X/PROD
1074 PROD = PROD-2.0D0
1075 140 CONTINUE
1076 160 ANSWR = FACT*SUM
1077 200 RETURN
1078 END
1079 C Seite 23
1080
1081
1082 SUBROUTINE BN(PCKR,NRINK,BSSLSP,CNEUMN)
1083 DIMENSION BSSLSP(NRINK),CNEUMN(NRINK)
1084 DOUBLE PRECISION PCKR,ANSWR,ANSA,ANSB,ANSC,CONN,CMULN,SNSA,SNSB,
1085 1SNSC,BSSLSP,CNEUMN
1086 NRANKI = NRINK
1087 NRANK = NRANKI-1
1088 NVAL = NRANK-1
1089 DO 40 I = 1,4
1090 CALL BESSEL(NVAL,PCKR,ANSWR,IERROR)
1091 IF(IERROR)20,20,32
1092 20 ANSA = ANSWR
1093 NVAL = NVAL+1
1094 CALL BESSEL(NVAL,PCKR,ANSWR,IERROR)
1095 IF(IERROR)24,24,28
1096 24 ANSB = ANSWR
1097 GO TO 60
1098 28 NVAL = NVAL-1
1099 32 NVAL = NVAL+NRANK
1100 40 CONTINUE
1101 STOP
1102 60 IF(NVAL-NRANK)100,100,64
1103 64 IEND = NVAL-NRANK
1104 CONN = DFLOAT(2*(NVAL-1)+1)
1105 DO 72 IP = 1,IEND
1106 ANSC = CONN*ANSA/PCKR-ANSB
1107 CONN = CONN-2.0D0
1108 ANSB = ANSA
1109 ANSA = ANSC
1110 72 CONTINUE
1111 100 BSSLSP(NRANKI) = ANSB
1112 BSSLSP(NRANKI-1) = ANSA
1113 CONN = DFLOAT(NRANK+NRANK-1)
1114 IE = NRANKI-2
1115 JE = IE
1116 DO 120 JB = 1,JE
1117 ANSC = CONN*ANSA/PCKR-ANSB
1118 BSSLSP(IE) = ANSC
1119 ANSB = ANSA
1120 ANSA = ANSC
1121 IE = IE-1
1122 CONN = CONN-2.0D0
1123 120 CONTINUE
1124 CMULN = 3.0D0
1125 SNSA = -DCOS(PCKR)/PCKR
1126 SNSB = -DCOS(PCKR)/(PCKR*PCKR)-DSIN(PCKR)/PCKR
1127 CNEUMN(1) = SNSA
1128 CNEUMN(2) = SNSB
1129 DO 280 I = 3,NRANKI
1130 SNSC = CMULN*SNSB/PCKR-SNSA
1131 CNEUMN(I) = SNSC
1132 SNSA = SNSB
1133 SNSB = SNSC
1134 CMULN = CMULN+2.0D0
1135 280 CONTINUE
1136 RETURN
1137 END
1138 C Seite 24
1139
1140
1141 SUBROUTINE CBESS (NORDER,ARGMNT,ANSWR,IERROR)
1142 COMPLEX*16 ARGMNT,ANSWR,X,SUM,APR,TOPR,CI,CNI,ACR,PROD,FACT
1143 IERROR = 0
1144 N = NORDER
1145 X = ARGMNT
1146 SUM = (1.0D0,0.0D0)
1147 APR = (1.0D0,0.0D0)
1148 TOPR = -(0.5D0,0.0D0)*X*X
1149 CI = (1.0D0,0.0D0)
1150 CNI = DCMPLX(DFLOAT(2*N+3),0.0D0)
1151 DO 60 I = 1,100
1152 ACR = TOPR*APR/(CI*CNI)
1153 SUM = SUM+ACR
1154 IF(CDABS(ACR/SUM)-1.0D-20)100,100,40
1155 40 APR = ACR
1156 CI = CI+(1.0D0,0.0D0)
1157 CNI = CNI+(2.0D0,0.0D0)
1158 60 CONTINUE
1159 IERROR = 1
1160 PRINT 10
1161 10 FORMAT('0 ERROR IN SUM OF CBESS')
1162 GO TO 200
1163 100 PROD = DCMPLX(DFLOAT(2*N+1),0.0D0)
1164 FACT = (1.0D0,0.0D0)
1165 IF(N)160,160,120
1166 120 DO 140 IFCT = 1,N
1167 FACT = FACT*X/PROD
1168 PROD = PROD-(2.0D0,0.0D0)
1169 140 CONTINUE
1170 160 ANSWR = FACT*SUM
1171 200 RETURN
1172 END
1173 C Seite 25
1174
1175
1176 SUBROUTINE CBN(PCKR,NRINK,BSSLSP,CNEUMN)
1177 DIMENSION BSSLSP(NRINK),CNEUMN(NRINK)
1178 COMPLEX*16 PCKR,ANSWR,ANSA,ANSB,ANSC,CONN,CMULN,SNSA,SNSB,SNSC,
1179 1BSSLSP,CNEUMN
1180 NRANKI = NRINK
1181 NRANK = NRANKI-1
1182 NVAL = NRANK-1
1183 DO 40 I = 1,4
1184 CALL CBESS(NVAL,PCKR,ANSWR,IERROR)
1185 IF(IERROR)20,20,32
1186 20 ANSA = ANSWR
1187 NVAL = NVAL+1
1188 CALL CBESS(NVAL,PCKR,ANSWR,IERROR)
1189 IF(IERROR)24,24,28
1190 24 ANSB = ANSWR
1191 GO TO 60
1192 28 NVAL = NVAL-1
1193 32 NVAL = NVAL+NRANK
1194 40 CONTINUE
1195 STOP
1196 60 IF(NVAL-NRANK)100,100,64
1197 64 IEND = NVAL-NRANK
1198 CONN = DCMPLX(DFLOAT(2*(NVAL-1)+1),0.0D0)
1199 DO 72 IP = 1,IEND
1200 ANSC = CONN*ANSA/PCKR-ANSB
1201 CONN = CONN-(2.0D0,0.0D0)
1202 ANSB = ANSA
1203 ANSA = ANSC
1204 72 CONTINUE
1205 100 BSSLSP(NRANKI) = ANSB
1206 BSSLSP(NRANKI-1) = ANSA
1207 CONN = DCMPLX(DFLOAT(NRANK+NRANK-1),0.0D0)
1208 IE = NRANKI-2
1209 JE = IE
1210 DO 120 JB = 1,JE
1211 ANSC = CONN*ANSA/PCKR-ANSB
1212 BSSLSP(IE) = ANSC
1213 ANSB = ANSA
1214 ANSA = ANSC
1215 IE = IE-1
1216 CONN = CONN-(2.0D0,0.0D0)
1217 120 CONTINUE
1218 CMULN = (3.0D0,0.0D0)
1219 SNSA = -CDCOS(PCKR)/PCKR
1220 SNSB = -CDCOS(PCKR)/(PCKR*PCKR)-CDSIN(PCKR)/PCKR
1221 CNEUMN(1) = SNSA
1222 CNEUMN(2) = SNSB
1223 DO 280 I = 3,NRANKI
1224 SNSC = CMULN*SNSB/PCKR-SNSA
1225 CNEUMN(I) = SNSC
1226 SNSA = SNSB
1227 SNSB = SNSC
1228 CMULN = CMULN+(2.0D0,0.0D0)
1229 280 CONTINUE
1230 RETURN
1231 END
1232 C Seite 26
1233
1234
1235 SUBROUTINE LEG(THETA,M,NRJNK,PNMLLG)
1236 DIMENSION PNMLLG(NRJNK)
1237 DOUBLE PRECISION THETA,PLA,PLB,PLC,DTWM,CNM,CNMM,CNMUL,PRODM,
1238 1QUANM,PNMLLG
1239 NRANKI = NRJNK
1240 KMV = M
1241 DTWM = DFLOAT (2*M+1)
1242 IF(THETA)16,4,16
1243 4 IF(KMV)12,12,6
1244 6 DO 8 ILG = 1,NRANKI
1245 PNMLLG(ILG) = 0.0D0
1246 8 CONTINUE
1247 GO TO 88
1248 12 PNMLLG(1) = 1.0D0
1249 PLA = 1.0D0
1250 GO TO 48
1251 16 IF(KMV)20,20,40
1252 20 PLA = 1.0D0
1253 PLB = DCOS(THETA)*PLA
1254 PNMLLG(1) = PLA
1255 PNMLLG(2) = PLB
1256 IBEG = 3
1257 GO TO 60
1258 40 DO 44 ILG = 1,KMV
1259 PNMLLG(ILG) = 0.0D0
1260 44 CONTINUE
1261 PRODM = 1.0D0
1262 QUANM = DFLOAT(KMV)
1263 DO 52 IFCT = 1,KMV
1264 QUANM = QUANM+1.0D0
1265 PRODM = QUANM*PRODM/2.0D0
1266 52 CONTINUE
1267 PLA = PRODM*DSIN(THETA)**KMV
1268 PNMLLG(KMV+1) = PLA
1269 48 PLB = DTWM*DCOS(THETA)*PLA
1270 PNMLLG(KMV+2) = PLB
1271 IBEG = KMV+3
1272 60 CNMUL = DFLOAT(IBEG+IBEG-3)
1273 CNM = 2.0D0
1274 CNMM = DTWM
1275 DO 80 ILGR = IBEG,NRANKI
1276 PLC = (CNMUL*DCOS(THETA)*PLB-CNMM*PLA)/CNM
1277 PNMLLG(ILGR) = PLC
1278 PLA = PLB
1279 PLB = PLC
1280 CNMUL = CNMUL+2.0D0
1281 CNM = CNM+1.0D0
1282 CNMM = CNMM+1.0D0
1283 80 CONTINUE
1284 88 RETURN
1285 END
1286 C Seite 27
1287
1288
1289 FUNCTION TRIXJ(J1,J2,J,M,FCT,N,L)
1290 IMPLICIT REAL*8(A-H,O-Z)
1291 DIMENSION FCT(1)
1292 INTEGER Z,ZMIN,ZMAX
1293 Y=1.0D0
1294 CC=0.0D0
1295 JSUM=J1+J2+J
1296 JM1=J1-IABS(M)
1297 JM2=J2-IABS(M)
1298 IF((MOD(JSUM,2).NE.0).OR.(MOD(JM1,2).NE.0).OR.(MOD(JM2,2).NE.0)
1299 1.OR.(JM1.LT.0).OR.(JM2.LT.0))
1300 2GO TO 3
1301 IF((J.GT.J1+J2).OR.(L.LT.IABS(J1-J2))) GO TO 1
1302 ZMIN=0
1303 IF(J-J2+M.LT.0) ZMIN=-J+J2-M
1304 IF(J-J1+M+ZMIN.LT.0) ZMIN=-J+J1-M
1305 ZMAX=J1+J2-J
1306 IF(J2-M-ZMAX.LT.0) ZMAX=J2-M
1307 IF(J1-M-ZMAX.LT.0) ZMAX=J1-M
1308 JA=(J1+M)/2+1
1309 JB=JA-M
1310 JC=(J2-M)/2+1
1311 JD=JC+M
1312 JE=J/2+1
1313 JF=(J1+J2-J)/2+1
1314 JG=JA+JB-JF
1315 JH=JC+JD-JF
1316 JJ=2*JE+JF-1
1317 IF(JJ.GT.N) Y=0.1D0**L
1318 IF(FCT(JJ))7,5,7
1319 7 CONTINUE
1320 IA=ZMIN/2
1321 IB=JF-IA+1
1322 IC=JB-IA+1
1323 ID=JC-IA+1
1324 IE=JA-JF+IA
1325 IF=JD-JF+IA
1326 FASE=1.0D0
1327 IF(MOD(IA,2).EQ.0) FASE=-FASE
1328 Z=ZMIN
1329 2 IA=IA+1
1330 IB=IB-1
1331 IC=IC-1
1332 ID=ID-1
1333 IE=IE+1
1334 IF=IF+1
1335 FASE=-FASE
1336 CC=CC+FASE/(FCT(IA)*FCT(IB)*FCT(IC)*FCT(ID)*FCT(IE)*FCT(IF))
1337 Z=Z+2
1338 IF(Z.LE.ZMAX) GO TO 2
1339 FASE=-DSIGN(1.0D0,CC)
1340 IF(MOD(J1-J2,4).EQ.0) FASE=-FASE
1341 CC=FASE*DSQRT(Y*FCT(JB)*FCT(JC)*FCT(JE)*CC*FCT(JA)*FCT(JD)*FCT(JE)
1342 1*CC*FCT(JF)*FCT(JG)*FCT(JH)/FCT(JJ))
1343 1 TRIXJ=CC
1344 RETURN
1345 3 PRINT 4
1346 4 FORMAT('0 ERROR IN ARGUMENT OF TRIXJ')
1347 CALL EXIT
1348 5 PRINT 6
1349
1350
1351 6 FORMAT('0 ERROR FACTORIALS EXCEEDED IN TIXJ')
1352 CALL EXIT
1353 END
1354 C Seite 28
1355
1356
1357 SUBROUTINE VPSI(RAD,THETA,FHI,NP,NB,M,ND,PSIR,PSITH,PSIFH,PN,U,V)
1358 DIMENSION PSIR(1),PSITH(1),PSIFH(1),PN(1),U(1),V(1)
1359 DOUBLE PRECISION RAD,R,THETA,T,FHI,F,PN,U,V,FCT,FR,FR1,FR2,B,C
1360 COMPLEX*16 FC,FC1,FC2,PSIR,PSITH,PSIFH
1361 COMMON/FAC/FCT(100),NRISK,NTRIX
1362 R = RAD
1363 T = THETA
1364 F = FHI
1365 NP1 = NP-1
1366 K = M
1367 N = ND
1368 L = N/2
1369 L1 = L+1
1370 L2 = L+2
1371 FC = (0.0D0,1.0D0)
1372 IF(NB.EQ.1) GO TO 5
1373 CALL BN(R,L1,U,V)
1374 B = DABS(R*R*(U(2)*V(1)-U(1)*V(2))-1.0D0)
1375 C = DABS(R*R*(U(L1)*V(L)-U(L)*V(L1))-1.0D0)
1376 PRINT 3
1377 3 FORMAT('0 BESSEL- NEUMAN- TEST FOR VPSI')
1378 PRINT 4,B,C
1379 4 FORMAT(2D25.15)
1380 5 CONTINUE
1381 CALL LEG(T,K,L2,PN)
1382 DO 42 I = 1,L
1383 I1 = I+1
1384 I2 = I+2
1385 IF(I-K)6,7,7
1386 6 FR = 0.0D0
1387 GO TO 10
1388 7 IF(K)8,8,9
1389 8 FR = DSQRT(DFLOAT(2*I+1)/(DFLOAT(I*I1)*16.0D0*DATAN(1.0D0)))
1390 GO TO 10
1391 9 FR = DSQRT(DFLOAT(2*I+1)*FCT(I1-K)/(DFLOAT(I*I1)*FCT(I1+K)*8.0D0*
1392 1DATAN(1.0D0)))
1393 10 CONTINUE
1394 IF(T)16,16,19
1395 16 CONTINUE
1396 IF(K-1)18,17,18
1397 17 CONTINUE
1398 FR1 = DSQRT(DFLOAT(2*I+1)/(32.0D0*DATAN(1.0D0)))
1399 FR2 = DSQRT(DFLOAT(2*I+1)/(32.0D0*DATAN(1.0D0)))
1400 GO TO 20
1401 18 CONTINUE
1402 FR1 = 0.0D0
1403 FR2 = 0.0D0
1404 GO TO 20
1405 19 CONTINUE
1406 FR1 = FR*PN(I1)*DFLOAT(K)/DSIN(T)
1407 FR2 = -FR*(DFLOAT(I1)*DCOS(T)*PN(I1)-DFLOAT(I1-K)*PN(I2))/DSIN(T)
1408 20 CONTINUE
1409 GO TO (30,31,32),NB
1410 30 CONTINUE
1411 FC1 = (-FC)**I1
1412 FC2 = (-FC)**I
1413 GO TO 33
1414 31 CONTINUE
1415 C Seite 29
1416
1417
1418 FC1 = DCMPLX(U(I1),0.0D0)
1419 FC2 = DCMPLX(U(I)-DFLOAT(I)*U(I1)/R,0.0D0)
1420
1421
1422 GO TO 33
1423 32 CONTINUE
1424 FC1 = DCMPLX(U(I1),V(I1))
1425 FC2 = DCMPLX(U(I)-DFLOAT(I)*U(I1)/R,V(I)-DFLOAT(I)*V(I1)/R)
1426 33 CONTINUE
1427 PSIR(2*I-1) = (0.0D0,0.0D0)
1428 IF(NB-1)34,34,35
1429 34 CONTINUE
1430 PSIR(2*I) = (0.0D0,0.0D0)
1431 35 CONTINUE
1432 GO TO (36,39),NP1
1433 36 CONTINUE
1434 IF(NB-1)38,38,37
1435 37 CONTINUE
1436 PSIR(2*I) = FC1*DCMPLX(DFLOAT(I*I1)*FR*PN(I1)*DSIN(DFLOAT(K)*F)/R,
1437 10.0D0)
1438 38 CONTINUE
1439 PSITH(2*I-1) = -FC1*DCMPLX(FR1*DSIN(DFLOAT(K)*F),0.0D0)
1440 PSITH(2*I) = FC2*DCMPLX(FR2*DSIN(DFLOAT(K)*F),0.0D0)
1441 PSIFH(2*I-1) = -FC1*DCMPLX(FR2*DCOS(DFLOAT(K)*F),0.0D0)
1442 PSIFH(2*I) = FC2*DCMPLX(FR1*DCOS(DFLOAT(K)*F),0.0D0)
1443 GO TO 42
1444 39 CONTINUE
1445 IF(NB-1)41,41,40
1446 40 CONTINUE
1447 PSIR(2*I) = FC1*DCMPLX(DFLOAT(I*I1)*FR*PN(I1)*DCOS(DFLOAT(K)*F)/R,
1448 10.0D0)
1449 41 CONTINUE
1450 PSITH(2*I-1) = FC1*DCMPLX(FR1*DCOS(DFLOAT(K)*F),0.0D0)
1451 PSITH(2*I) = FC2*DCMPLX(FR2*DCOS(DFLOAT(K)*F),0.0D0)
1452 PSIFH(2*I-1) = -FC1*DCMPLX(FR2*DSIN(DFLOAT(K)*F),0.0D0)
1453 PSIFH(2*I) = -FC2*DCMPLX(FR1*DSIN(DFLOAT(K)*F),0.0D0)
1454 42 CONTINUE
1455 RETURN
1456 END
1457 C Seite 30
1458
1459
1460 SUBROUTINE VKOEF (BHETA,NP,M,ND,AP,PN)
1461 DIMENSION PN(1),AP(1)
1462 DOUBLE PRECISION BHETA,U,DI,PN,FR,FCT
1463 COMPLEX*16 FC1,FC2,fc3,AP
1464 COMMON/FAC/FCT(100),NRISK,NTRIX
1465 N = NP
1466 N1 = N+1
1467 K = M
1468 L = ND/2
1469 L2 = L+2
1470 U = BHETA
1471 DI = DATAN(1.0D0)
1472 FC1 = (0.0D0,1.0D0)
1473 CALL LEG(U,K,L2,PN)
1474 IF(U)10,10,20
1475 10 DO 100 I = 1,L
1476 IF(K-1)11,12,11
1477 11 AP(2*I-1) = (0.0D0,0.0D0)
1478 AP(2*I) = (0.0D0,0.0D0)
1479 GO TO 100
1480 12 FR = DSQRT(8.0D0*DI*DFLOAT(2*I+1))
1481 FC2 = DCMPLX(FR,0.0D0)
1482 GO TO (13,14),N1
1483 13 AP(2*I-1) = -(FC1**I)*FC2
1484 AP(2*I) = -(FC1**(I+1))*FC2
1485 GO TO 100
1486 14 AP(2*I-1) = (FC1**I)*FC2
1487 AP(2*I) = (FC1**(I+3))*FC2
1488 100 CONTINUE
1489 GO TO 300
1490 20 DO 200 I = 1,L
1491 I1 = I+1
1492 I2 = I+2
1493 IF(K)21,21,30
1494 21 FR = 4.0D0*DSQRT((DI*DFLOAT((2*I+1)*(I+1)))/DFLOAT(I))*(DCOS(U)*
1495 1PN(I1)-PN(I2))/DSIN(U)
1496 FC2 = DCMPLX(FR,0.0D0)
1497 GO TO (22,23),N1
1498 22 AP(2*I-1) = (FC1**I)*FC2
1499 AP(2*I) = (0.0D0,0.0D0)
1500 GO TO 200
1501 23 AP(2*I-1) = (0.0D0,0.0D0)
1502 AP(2*I) = (FC1**(I+1))*FC2
1503 GO TO 200
1504 30 IF(I-K)34,31,31
1505 31 FR = 4.0D0*DSQRT((2.0D0*DI*DFLOAT(2*I+1)*FCT(I-K+1))/(DFLOAT(I*(I
1506 1+1))*FCT(I+K+1)))
1507 FC2 = DCMPLX((FR*(DFLOAT(I+1)*DCOS(U)*PN(I1)-DFLOAT(I-K+1)*PN(I2)
1508 1))/DSIN(U),0.0D0)
1509 FC3 = DCMPLX((DFLOAT(K)*FR*PN(I1))/DSIN(U),0.0D0)
1510 GO TO (32,33),N1
1511 32 AP(2*I-1) = (FC1**I)*FC2
1512 AP(2*I) = -(FC1**(I+1))*FC3
1513 GO TO 200
1514 33 AP(2*I-1) = (FC1**I)*FC3
1515 AP(2*I) = (FC1**(I+1))*FC2
1516 GO TO 200
1517 34 AP(2*I-1) = (0.0D0,0.0D0)
1518 AP(2*I) = (0.0D0,0.0D0)
1519 200 CONTINUE
1520
1521
1522 300 RETURN
1523 END
1524 C Seite 31
1525
1526
1527 SUBROUTINE VR(AR,NP,M,ND,R1,R2,R3,R4,X1,X2,Y1,Y2)
1528 DIMENSION R1(ND,1),R2(ND,1),R3(ND,1),R4(ND,1),X1(1),X2(1),Y1(1),Y2
1529 1(1)
1530 DOUBLE PRECISION AR,U,V,F3J1,F3J2,F3J3,FACT1,FACT2,FACT3,FACT4,
1531 1FACT5,FACT6,FACT7,FACT8,FACT9,FACT10,ACR1,ACR2,ACR3,ACR4,X1,X2,
1532 1Y1,Y2,B1,B2,CN1,CN2,FCT
1533 COMPLEX*16 R1,R2,R3,R4,W1,W2,W3,W4,W5,W6,W7,W8
1534 COMMON/FAC/FCT(100),NRISK,NTRIX
1535 NP1 = NP+1
1536 N = M
1537 NK = ND+1
1538 U = AR
1539 V = 0.5D0*AR
1540 CALL BN(U,NK,X1,Y1)
1541 CALL BN(V,NK,X2,Y2)
1542 L1 = NK
1543 L = NK-1
1544 B1 = DABS(U**2*(X1(2)*Y1(1)-X1(1)*Y1(2))-1.0D0)
1545 B2 = DABS(V**2*(X1(2)*Y2(1)-X2(1)*Y2(2))-1.0D0)
1546 CN1 = DABS(U**2*(X1(L1)*Y1(L)-X1(L)*Y1(L1))-1.0D0)
1547 CN2 = DABS(V**2*(X2(L1)*Y2(L)-X2(L)*Y2(L1))-1.0D0)
1548 PRINT 89
1549 89 FORMAT('0 BESSEL- NEUMAN- TEST FOR VR')
1550 PRINT 90, B1,CN1,B2,CN2
1551 90 FORMAT(4D25.15)
1552 DO 7 I = 1,L
1553 DO 7 J = 1,L
1554 R1(I,J) = (0.0D0,0.0D0)
1555 R2(I,J) = (0.0D0,0.0D0)
1556 R3(I,J) = (0.0D0,0.0D0)
1557 R4(I,J) = (0.0D0,0.0D0)
1558 7 CONTINUE
1559 NR = L/2
1560 DO 200 I = 1,NR
1561 DO 200 J = 1,NR
1562 IF(N-I)8,8,200
1563 8 IF(N-J)9,9,200
1564 9 L1 = I+J+1
1565 W1 = (0.0D0,0.0D0)
1566 W2 = (0.0D0,0.0D0)
1567 W3 = (0.0D0,0.0D0)
1568 W4 = (0.0D0,0.0D0)
1569 W5 = (0.0D0,0.0D0)
1570 W6 = (0.0D0,0.0D0)
1571 W7 = (0.0D0,0.0D0)
1572 W8 = (0.0D0,0.0D0)
1573 DO 100 L = 1,L1
1574 K = L-1
1575 IF(K-IABS(I-J))100,10,10
1576 10 CONTINUE
1577 IF(N)11,11,12
1578 11 IF(MOD((I+J+K),2).NE.0) GO TO 100
1579 FACT1 = 1.0D0
1580 GO TO 13
1581 12 FACT1 = DFLOAT((-1)**N)
1582 13 J1 = 2*I
1583 J2 = 2*J
1584 J3 = 2*K
1585 M1 = 2*N
1586 C Seite 32
1587
1588
1589 F3J1 = TRIXJ(J1,J2,J3,M1,FCT,NRISK,NTRIX)
1590 FACT2 = 0.5D0*DFLOAT(2*K+1)*F3J1*
1591
1592
1593 1DSQRT(DFLOAT((2*I+1)*(2*J+1))/DFLOAT(I*(I+1)*J*(J+1)))
1594 FACT3 = FACT1*FACT2
1595 IF(MOD((J-I+K),2).NE.0) GO TO 27
1596 IF(J-I+K)25,23,24
1597 23 FACT4 = 1.0D0
1598 GO TO 26
1599 24 FACT4 = DFLOAT((-1)**((J-I+K)/2))
1600 GO TO 26
1601 25 FACT4 = DFLOAT((-1)**((I-J-K)/2))
1602 26 J1 = 2*I
1603 J2 = 2*J
1604 J3 = 2*K
1605 M1 = 0
1606 F3J2 = TRIXJ(J1,J2,J3,M1,FCT,NRISK,NTRIX)
1607 FACT5 = DFLOAT(I*(I+1)+J*(J+1)-K*(K+1))*F3J2
1608 GO TO 28
1609 27 FACT6 = 0.0D0
1610 GO TO 29
1611 28 FACT6 = FACT4*FACT5
1612 29 IF(K-IABS(I-J)-1)44,30,30
1613 30 IF(MOD((J-I+K+1),2).NE.0) GO TO 44
1614 IF(J-I+K+1)33,31,32
1615 31 FACT7 = 1.0D0
1616 GO TO 41
1617 32 FACT7 = DFLOAT((-1)**((J-I+K+1)/2))
1618 GO TO 41
1619 33 FACT7 = DFLOAT((-1)**((I-J-K-1)/2))
1620 41 J1 = 2*I
1621 J2 = 2*J
1622 J3 = 2*(K-1)
1623 M1 = 0
1624 F3J3 = TRIXJ(J1,J2,J3,M1,FCT,NRSIK,NTRIX)
1625 IF(IABS(I-J))42,42,43
1626 42 FACT8 = DFLOAT(K)*DSQRT(DFLOAT((I+J-1)**2-K**2))*F3J3
1627 GO TO 45
1628 43 FACT8 = DSQRT(DFLOAT((K**2-IABS(I-J)**2)*((I+J+1)**2-K**2)))*F3J3
1629 GO TO 45
1630 44 FACT9 = 0.0D0
1631 GO TO 50
1632 45 FACT9 = FACT7*FACT8
1633 50 CONTINUE
1634 IF(K)51,51,52
1635 51 FACT10 = 1.0D0
1636 GO TO 53
1637 52 FACT10 = DFLOAT((-1)**K)
1638 53 ACR1 = FACT3*FACT6
1639 ACR2 = FACT10*ACR1
1640 ACR3 = FACT3*FACT9
1641 ACR4 = FACT10*ACR3
1642 K1 = K+1
1643 W1 = W1+DCMPLX(ACR1*X1(K1),0.0D0)
1644 W2 = W2+DCMPLX(ACR1*X1(K1),ACR1*Y1(K1))
1645 W3 = W3+DCMPLX(ACR2*X1(K1),ACR2*Y1(K1))
1646 W4 = W4+DCMPLX(ACR1*X2(K1),0.0D0)
1647 IF(N)100,100,99
1648 C Seite 33
1649
1650
1651 99 CONTINUE
1652 W5 = W5+DCMPLX(ACR3*X1(K1),0.0D0)
1653 W6 = W6+DCMPLX(ACR3*X1(K1),ACR3*Y1(K1))
1654 W7 = W7+DCMPLX(ACR4*X1(K1),ACR4*Y1(K1))
1655
1656
1657 W8 = W8+DCMPLX(ACR3*X2(K1),0.0D0)
1658 100 CONTINUE
1659 R1(2*I-1,2*J-1) = W1
1660 R2(2*I-1,2*J-1) = W2
1661 R3(2*I-1,2*J-1) = W3
1662 R4(2*I-1,2*J-1) = W4
1663 R1(2*I,2*J) = W1
1664 R2(2*I,2*J) = W2
1665 R3(2*I,2*J) = W3
1666 R4(2*I,2*J) = W4
1667 GO TO (101,103),NP1
1668 101 CONTINUE
1669 IF(N)200,200,102
1670 102 CONTINUE
1671 R1(2*I-1,2*J) = W5
1672 R2(2*I-1,2*J) = W6
1673 R3(2*I-1,2*J) = W7
1674 R4(2*I-1,2*J) = W8
1675 R1(2*I,2*J-1) = -W5
1676 R2(2*I,2*J-1) = -W6
1677 R3(2*I,2*J-1) = -W7
1678 R4(2*I,2*J-1) = -W8
1679 GO TO 200
1680 103 CONTINUE
1681 IF(N)200,200,104
1682 104 CONTINUE
1683 R1(2*I-1,2*J) = -W5
1684 R2(2*I-1,2*J) = -W6
1685 R3(2*I-1,2*J) = -W7
1686 R4(2*I-1,2*J) = -W8
1687 R1(2*I,2*J-1) = W5
1688 R2(2*I,2*J-1) = W6
1689 R3(2*I,2*J-1) = W7
1690 R4(2*I,2*J-1) = W5
1691 200 CONTINUE
1692 RETURN
1693 END
1694 C Seite 34
1695
1696
1697 SUBROUTINE MCNV(A,N,D,L,M)
1698 DIMENSION A(1),L(1),M(1)
1699 COMPLEX*16 A,D,BIGA,HOLD
1700 D = (1.0D0,0.0D0)
1701 NK=-N
1702 DO 80 K=1,N
1703 NK=NK+N
1704 L(K)=K
1705 M(K)=K
1706 KK=NK+K
1707 BIGA=A(KK)
1708 DO 20 J=K,N
1709 IZ=N*(J-1)
1710 DO 20 I=K,N
1711 IJ=IZ+1
1712 10 IF(CDABS(BIGA)-CDABS(A(IJ))) 15,20,20
1713 15 BIGA=A(IJ)
1714 L(K)=I
1715 M(K)=J
1716 20 CONTINUE
1717 J=L(K)
1718 IF(J-K) 35,35,25
1719 25 KI=K-N
1720 DO 30 I=1,N
1721 KI=KI+N
1722 HOLD=-A(KI)
1723 JI=KI-K+J
1724 A(KI)=A(JI)
1725 30 A(JI) =HOLD
1726 35 I=M(K)
1727 IF(I-K) 45,45,38
1728 38 JP=N*(I-1)
1729 DO 40 J=1,N
1730 JK=NK+J
1731 JI=JP+J
1732 HOLD=-A(JK)
1733 A(JK)=A(JI)
1734 40 A(JI) =HOLD
1735 45 IF(CDABS(BIGA)) 48,46,48
1736 46 D = (0.0D0,0.0D0)
1737 RETURN
1738 48 DO 55 I=1,N
1739 IF(I-K) 50,55,50
1740 50 IK=NK+1
1741 A(IK)=A(IK)/(-BIGA)
1742 55 CONTINUE
1743 DO 65 I=1,N
1744 IK=NK+1
1745 HOLD=A(IK)
1746 IJ=I-N
1747 DO 65 J=1,N
1748 IJ=IJ+N
1749 IF(I-K) 60,65,60
1750 60 IF(J-K) 62,65,62
1751 62 KJ=IJ-I+K
1752 A(IJ)=HOLD*A(KJ)+A(IJ)
1753 65 CONTINUE
1754 KJ=K-N
1755 DO 75 J=1,N
1756 C Seite 35
1757
1758
1759 KJ=KJ+N
1760
1761
1762 IF(J-K) 70,75,70
1763 70 A(KJ)=A(KJ)/BIGA
1764 75 CONTINUE
1765 D=D*BIGA
1766 A(KK) = (1.0D0,0.0D0)/BIGA
1767 80 CONTINUE
1768 K=N
1769 100 K=(K-1)
1770 IF(K) 150,150,105
1771 105 I=L(K)
1772
1773
1774 IF(I-K) 120,120,108
1775 108 JQ=N*(K-1)
1776 JR=N*(I-1)
1777 DO 110 J=1,N
1778 JK=JQ+J
1779 HOLD=A(JK)
1780 JI=JR+J
1781 A(JK)=-A(JI)
1782 110 A(JI) =HOLD
1783 120 J=M(K)
1784 IF(J-K) 100,100,125
1785 125 KI=K-N
1786 DO 130 I=1,N
1787 KI=KI+N
1788 HOLD=A(KI)
1789 JI=KI-K+J
1790 A(KI)=-A(JI)
1791 130 A(JI) =HOLD
1792 GO TO 100
1793 150 RETURN
1794 END
1795 C Seite 36
1796
1797
1798 SUBROUTINE COND(M,ND,RE,IM)
1799 DIMENSION RE(ND,1),IM(ND,1)
1800 DOUBLE PRECISION RE,IM,SCALE
1801 MD = 1
1802 IF(M.GT.1) MD = 2*M-1
1803 MD1 = MD+1
1804 NBGR = ND
1805 NROW = NBGR
1806 DO 60 KR = MD1,NBGR
1807 SCALE = 1.0D0/IM(NROW,NROW)
1808 DO 8 LC = MD,NBGR
1809 RE(NROW,LC) = SCALE*RE(NROW,LC)
1810 IM(NROW,LC) = SCALE*IM(NROW,LC)
1811 8 CONTINUE
1812 MROW = NROW-1
1813 DO 20 MR = MD,MROW
1814 SCALE = IM(MR,NROW)
1815 DO 16 MC = MD,NBGR
1816 RE(MR,MC) = RE(MR,MC)-SCALE*RE(NROW,MC)
1817 IM(MR,MC) = IM(MR,MC)-SCALE*IM(NROW,MC)
1818 16 CONTINUE
1819 20 CONTINUE
1820 NROW = NROW-1
1821 60 CONTINUE
1822 NROW = NBGR-1
1823 DO 80 I = 1,NROW
1824 IB = I+1
1825 DO 72 J = IB,NBGR
1826 IM(I,J) = 0.0D0
1827 72 CONTINUE
1828 80 CONTINUE
1829 RETURN
1830 END
1831 C Seite 37
1832
1833
1834 SUBROUTINE ORTHO(M,ND,RE,IM,X,Y)
1835 DIMENSION RE(ND,1),IM(ND,1),X(1),Y(1)
1836 DOUBLE PRECISION RE,IM,X,Y,SUM1,SUM2
1837 MD = 1
1838 IF(M.GT.1) MD = 2*M-1
1839 NBGR = ND
1840 SUM1 = 0.0D0
1841 DO 20 K = MD,NBGR
1842 SUM1 = SUM1+RE(NBGR,K)**2+IM(NBGR,K)**2
1843 20 CONTINUE
1844 SUM1 = DSQRT(SUM1)
1845 DO 28 K = MD,NBGR
1846 RE(NBGR,K) = RE(NBGR,K)/SUM1
1847 IM(NBGR,K) = IM(NBGR,K)/SUM1
1848 28 CONTINUE
1849 NMI = NBGR-1
1850 NROW = NBGR
1851 DO 100 I = MD,NMI
1852 NROW = NROW-1
1853 MROW = NROW
1854 DO 36 K = MD,NBGR
1855 X(K) = RE(NROW,K)
1856 Y(K) = IM(NROW,K)
1857 36 CONTINUE
1858 DO 80 J = NROW,NMI
1859 SUM1 = 0.0D0
1860 SUM2 = 0.0D0
1861 MROW = MROW+1
1862 DO 40 K = MD,NBGR
1863 SUM1 = SUM1+RE(MROW,K)*RE(NROW,K)+IM(MROW,K)*IM(NROW,K)
1864 SUM2 = SUM2+RE(MROW,K)*IM(NROW,K)-IM(MROW,K)*RE(NROW,K)
1865 40 CONTINUE
1866 DO 48 K = MD,NBGR
1867 X(K) = X(K)-SUM1*RE(MROW,K)+SUM2*IM(MROW,K)
1868 Y(K) = Y(K)-SUM1*IM(MROW,K)-SUM2*RE(MROW,K)
1869 48 CONTINUE
1870 80 CONTINUE
1871 SUM1 = 0.0D0
1872 DO 84 K = MD,NBGR
1873 SUM1 = SUM1+X(K)**2+Y(K)**2
1874 84 CONTINUE
1875 SUM1 = DSQRT(SUM1)
1876 DO 88 K = MD,NBGR
1877 RE(NROW,K) = X(K)/SUM1
1878 IM(NROW,K) = Y(K)/SUM1
1879 88 CONTINUE
1880 100 CONTINUE
1881 RETURN
1882 END
1883 C Seite 38
1884
1885
1886 SUBROUTINE PERFT(NP,M,NR,ND,AR,AI,BR,BI,T,RE,IM,X,Y)
1887 DIMENSION AR(NR,1),AI(NR,1),BR(NR,1),BI(NR,1),T(ND,1),RE(ND,1),
1888 1IM(ND,1),X(1),Y(1)
1889 DOUBLE PRECISION AR,AI,BR,BI,RE,IM,X,Y,FAC
1890 COMPLEX*16 T
1891 MD = 1
1892 IF(M.GT.1) MD = 2*M-1
1893 NR = ND/2
1894 FAC = 1.0D0
1895 IF(NP.GT.0) FAC = -1.0D0
1896 DO 10 I = 1,NR
1897 DO 10 J = 1,NR
1898 RE(2*I-1,2*J-1) = AR(I,J)
1899 RE(2*I-1,2*J) = FAC*BR(I,J)
1900 RE(2*I,2*J-1) = FAC*BR(I,J)
1901 RE(2*I,2*J) = -AR(I,J)
1902 IM(2*I-1,2*J-1) = AI(I,J)
1903 IM(2*I-1,2*J) = FAC*BI(I,J)
1904 IM(2*I,2*J-1) = FAC*BI(I,J)
1905 IM(2*I,2*J) = -AI(I,J)
1906 IF(1.EQ.J) IM(2*I,2*I) = 1.0D0-AI(I,I)
1907 10 CONTINUE
1908 CALL COND(M,ND,RE,IM)
1909 CALL ORTHO(M,ND,RE,IM,X,Y)
1910 DO 11 I = 1,ND
1911 DO 11 J = 1,ND
1912 T(I,J) = (0.0D0,0.0D0)
1913 11 CONTINUE
1914 DO 12 I = MD,ND
1915 DO 12 J = MD,ND
1916 DO 12 K = MD,ND
1917 T(I,J) = T(I,J)-DCMPLX(RE(K,J)*RE(K,J),-IM(K,I)*RE(K,J))
1918 12 CONTINUE
1919 RETURN
1920 END
1921 C Seite 39
1922
1923
1924 SUBROUTINE LINE(THETA,NIP,C,ALPHA,R,DR)
1925 DOUBLE PRECISION THETA,C,ALPHA,R,DR
1926 R = C*DSIN(ALPHA)/DSIN(THETA+DFLOAT(NIP)*ALPHA)
1927 DR = -R/DTAN(THETA+DFLOAT(NIP)*ALPHA)
1928 RETURN
1929 END
1930
1931
1932 SUBROUTINE TRCIR(STH,CTH,C,A,R,DR)
1933 DOUBLE PRECISION STH,CTH,C,A,R,DR,X,ST,CT
1934 ST = C*STH
1935 CT = C*CTH
1936 X = DSQRT(A**2-ST**2)
1937 R = CT+X
1938 DR = -ST-ST*CT/X
1939 RETURN
1940 END
1941
1942
1943 SUBROUTINE ELLIPS(STH,CTH,AZ,BRA,R,DR)
1944 DOUBLE PRECISION STH,CTH,AZ,BRA,R,DR
1945 R = AZ/DSQRT (CTH**2+(AZ*STH/BRA)**2)
1946 DR = R**3*STH*CTH*(1.0D0-(AZ/BRA)**2)/AZ**2
1947 RETURN
1948 END
1949
1950
1951 SUBROUTINE TRELLI(STH,CTH,AZ,BRA,C,R,DR)
1952 DOUBLE PRECISION STH,CTH,AZ,BRA,C,R,DR,NE,RO
1953 NE = (AZ*STH)**2+(BRA*CTH)**2
1954 RO = NE-(C*STH)**2
1955 RO = DSQRT(RO)
1956 R = (BRA**2*C*CTH+AZ*BRA*RO)/NE
1957 DR = -2.0D0*STH*CTH*(AZ**2-BRA**2)*R/NE+(-BRA**2*C*STH+STH*CTH*AZ*
1958 1BRA*(AZ**2-BRA**2-C**2)/RO)/NE
1959 RETURN
1960 END
1961 C Seite 40
1962
1963 C Sphere-cone-sphere.
1964
1965 A = 1.0D0
1966 ALPHA = PI*15.0D0/180.0D0
1967 NDLTH(1) = 64
1968 NDLTH(2) = 64
1969 NDLTH(3) = 64
1970 CALL SPCOSP(ALPHA,A,THETA1,THETA2)
1971 NSECT = 3
1972 NIP = -1
1973 B = 0.5D0*A*(1.0D0-DSIN(ALPHA))/(1.0D0+DSIN(ALPHA))
1974 C = 0.5D0*A*(DSIN(ALPHA)**2+DSIN(ALPHA)+2.0D0)/(DSIN(ALPHA)*(1.0D0
1975 1+DSIN(ALPHA)))
1976 D = -0.5D0*A
1977 E = A/(1.0D0+DSIN(ALPHA))
1978 THETA3 = PI
1979 CDH(1) = THETA1/DFLOAT(NDLTH(1))
1980 CDH(2) = (THETA2-THETA1)/DFLOAT(NDLTH(2))
1981 CDH(3) = (THETA3-THETA2)/DFLOAT(NDLTH(3))
1982 GO TO (1,2,3),ISECT
1983 1 CONTINUE
1984 CALL TRCIR(STH,CTH,B,A,R,DR)
1985 GO TO 4
1986 2 CONTINUE
1987 CALL LINE(THETA,NIP,C,ALPHA,R,DR)
1988 GO TO 4
1989 3 CONTINUE
1990 CALL TRCIR(STH,CTH,D,E,R,DR)
1991 4 CONTINUE
1992
1993
1994 SUBROUTINE SPCOSP(ALPHA,A,THETA1,THETA2)
1995 DOUBLE PRECISION ALPHA,SNA,CSA,A,BDA,Q,THETA1,THETA2
1996 SNA = DSIN(ALPHA)
1997 CSA = DCOS(ALPHA)
1998 BDA = 1.0D0/(1.0D0+SNA)
1999 Q = (1.0D0-BDA)*(1.0D0-SNA)/2.0D0
2000 THETA1 = DATAN(SNA*CSA/(Q-SNA**2))
2001 THETA2 = DATAN(BDA*SNA*CSA/(1.0D0-Q-BDA*CSA**2))
2002 THETA2 = 4.0D0*DATAN(1.0D0)-THETA2
2003 RETURN
2004 END
2005 C Seite 41
2006
2007 C Translated circle
2008
2009 A = 1.0D0
2010 C = -0.1D0
2011 NDLTH(1) = 192
2012 NSECT = 1
2013 THETA1 = PI
2014 CDH(1) = THETA/DFLOAT(NDLTH(1))
2015
2016 CALL TRCIR(STH,CTH,C,A,R,DR)
2017
2018
2019
2020
2021
2022 C Ellips
2023
2024 A = 1.0D0
2025 B = 0.5D0
2026 NDLTH(1) = 1
2027 THETA1 = PI
2028 CDH(1) = THETA/DFLOAT(NDLTH(1))
2029
2030 CALL ELLIPS (STH,CTH,A,B,R,DR)
2031
2032
2033
2034 C Translated ellips
2035
2036 A = 1.0D0
2037 B = 0.5D0
2038 C = -0.1D0
2039 NDLTH(1) = 192
2040 NSECT = 1
2041 THETA1 = PI
2042 CDH(1) = THETA/DFLOAT(NDLTH(1))
2043
2044 CALL TRELLI(STH,CTH,A,B,C,R,DR)
2045 C Seite 42